#include <stdio.h>
#include <iostream>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include "eigen3/Eigen/IterativeLinearSolvers"
#include "eigen3/Eigen/Sparse"
using namespace Eigen;
int main()
{
  int dim,i,j,k,n;
  char c;
  double d,err;
  dim=119425;n=825737;
  clock_t start_time = clock();

  SparseMatrix<double> mat(dim,dim);
  mat.reserve(VectorXi::Constant(dim,30));

  FILE* f = fopen("spmat1.txt","r");    //read information
  for(k=0;k<n;k++)
    {
      fscanf(f,"%d%c%d%c%le",&i,&c,&j,&c,&d);
      mat.insert(i-1,j-1)=d;
    }
  fclose(f);
  
  mat.makeCompressed();
  VectorXd b(dim);
  f = fopen("spmat2.txt","r");
  for(k=0;k<dim;k++)
    {
      fscanf(f,"%le",&d);
      b(k)=d;
    }
  fclose(f);

  VectorXd x(dim);
  SparseLU<SparseMatrix<double>,COLAMDOrdering<int> >solver;    //solving with LU
  solver.analyzePattern(mat);
  solver.factorize(mat);
  x = solver.solve(b);
  VectorXd b1(dim);
  b1 = mat*x;
  VectorXd r(dim);
  r = b1-b;
  err = sqrt(r.dot(r));
  printf("the residual is : %le\n",err);
  clock_t  end_time = clock();
  printf("the running time is %lf\n",(double)(end_time - start_time)/CLOCKS_PER_SEC);
  return 0;
 
}
